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Abstract 

We investigate the dynamics of spinning binaries of compact objects at the next-to-leading 
order in the quadratic-in-spin effects, which corresponds to the third post-Newtonian order (3PN). 
Using a Dixon-type multipolar formalism for spinning point particles endowed with spin-induced 
quadrupoles and computing iteratively in harmonic coordinates the relevant pieces of the PN metric 
within the near zone, we derive the post-Newtonian equations of motion as well as the equations 
of spin precession. We find full equivalence with available results. We then focus on the far-zone 
field produced by those systems and obtain the previously unknown 3PN spin contributions to 
the gravitational-wave energy flux by means of the multipolar post-Minkowskian (MPM) wave 
generation formalism. Our results are presented in the center-of-mass frame for generic orbits, 
before being further specialized to the case of spin-aligned, circular orbits. We derive the orbital 
phase of the binary based on the energy balance equation and briefly discuss the relevance of the 
new terms. 

PACS numbers: 04.25.Nx, 04.25.dg, 04.30.-w, 97.80.-d, 97.60.Jd, 95.30.Sf 


* alejandro.bohe@aei.mpg.de 
f faye@iap.fr 
i smarsat@umd.edu 
§ porter@apc.univ-paris7.fr 


1 



I. INTRODUCTION 


Coalescing binary systems composed of stellar mass black-holes and/or neutron stars are 
among the most promising sources for a first direct detection of gravitational waves (GW) 
by the network of ground-based interferometers formed by GEO-HF [l) and the advanced 
version of the detectors LIGO [2j and Virgo |3j, which should resume their science runs from 
2015, approaching gradually their design sensitivity, expected to be better by an order of 
magnitude than that of the first generation. The cryogenic detector KAGRA ]U] will join 
them in a near future. Further ahead, the space-based observatory eLISA [ 5 ], 6j] — a serious 
proposal for the mission recently announced by the European Space Agency — will allow 
us to scan a different frequency band where we expect to detect, notably, GW emitted by 
supermassive black-hole binaries before merger. 

Extraction of the signal from the noisy data by means of matched filtering techniques 
and source parameter estimation both require an accurate modeling of the waveform. For 
binary systems of compact objects, the inspiralling phase of the coalescence can be modeled 
extremely well by resorting to the perturbative post-Newtonian (PN) scheme (see (?] for a 
review), in which all quantities of interest are expanded as formal series in powers of 1/c. 
For non-spinning (NS) systems, the phase of the waveform is currently known up to the 
order 3.5PN (i.e. including corrections up to 1/c 7 ), whereas the full polarizations have been 
obtained up to the order 3PN [§| (with the dominant quadrupole and octupole modes in the 
decomposition of the waveform in spin-weighted spherical harmonics known up to the order 

3.5PN@,0)- 

In recent years, motivated by astrophysical observations suggesting that black holes in 
our universe can have significant spins, considerable effort has been devoted to investigating 
higher order corrections to the spin effects in the binary dynamics, mostly restricted to 
the conservative piece of the body evolution in the near zone. While for the neutron stars 
observed so far, the largest dimensionless spin magnitude ever measured 0 is only y Q > 
(and may reasonably be assumed to be much smaller for typical expected observations), the 
spin of a black-hole might be commonly close to its maximal value [l2-T5||. Then, its effect 
on the waveform can be fairly strong and, in particular, for spins misaligned with the orbital 
angular momentum of the system, the dynamics becomes much more involved as the orbital 
plane undergoes precession, resulting in large modulations of the waveforms [l6|, 17]. Even 
in the simpler case where the spins are aligned with the orbital angular momentum, they 
significantly affect the inspiral rate of the binary, i.e. the frequency evolution of the signal, 
starting at the 1.5 PN order (see for instance Ref. [l8j for a detailed study of the effect of 
the spin on the waveform quantified in terms of figures of merit relevant to data analysis). 
To make all factors 1/c appear explicitly in this paper, we rescale the physical spin variable 

“^physical aS 

S c ^physical Gm y, (1 • 1) 


where y is the dimensionless spin, with value 1 for an extremal Kerr black hole. 

The calculation of the spin PN corrections to the conservative part of the dynamics and, 
to some extent, to the radiation field of the binary beyond the leading order contributions has 
been tackled using essentially three different approaches: (i) a Hamiltonian approach that 
strongly relies on the use of the (second) Arnowitt-Deser-Missner (ADM) gauge [l9j, and 
in which the dissipative part of the dynamics, demanding a special treatment, is generally 
discarded (see however Ref. (20[), (ii) an effective field theory (EFT) Lagrangian formal¬ 
ism j2ll, HH, whose application to binary systems in general relativity has been actively 
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developed since the mid-2000’s, and (in) a post-Newtonian iteration scheme in harmonic co¬ 
ordinates (PNISH), reviewed in Ref. |7|, which we follow in the present paper. The existence 
of those three independent methods permits important checks of calculations that are often 
tedious, whenever quantities are available at the same order in more than one formalism. 

The binary dynamics at the spin-orbit level (i.e. linear-in-spin effects, which will be 
referred to as SO from now on) are known up to the order 3.5PN in both the PNISH and 
ADM approaches [H-28|, and to the order 2.5PN in the EFT framework (29[ 30j. On the 
other hand, quadratic-in-spin corrections (labeled as SS throughout the paper) have been 


obtained to the order 2PN in the PNISH formalism 
EFT formalisms they are known up to the order 3PN 
,Sj S <2 interactions 
derived 



3311 . while in both the ADM and 
371 ]. and even 4PN for the simpler 
Higher-order-in-spin corrections have also been recently 


As for the spin contributions to the radiation held, they have mostly 


been computed by using the same usual combination of the MPM and PNISH approaches as 
in the present paper, although partial results required for the calculation of the 3PN flux 


and the 2.5PN waveform 44] have been obtained within the EFT approach. The en ergy 


flux of gravitational-wave radiation is known up to the order 4PN at the SO level 45 


47 


33 


whereas at the SS level only the leading order (2PN) terms were known until now 
Moreover, the leading order cubic-in-spin terms, which arise at 3.5PN, have been calculated 
very recently ] 42| . 

Our goal here will be to determine, within the PNISH approach, the 3PN (i.e. next-to- 
leading order) spin-spin corrections entering both the source dynamics (thereby providing 
an additional confirmation of the ADM and EFT results already available at this order) 
and most importantly the energy flux, thus completing the knowledge of all the spinning 
corrections to the phasing formula up to the 3PN order. At the next order 3.5PN, the 
only remaining unknown terms all come from a SS tail contribution. By contrast, the spin 
corrections to the full gravitational-wave polarizations are only known to the poorer 2PN 
and we postpone to future work the task of obtaining all the corrections 


33, 48 


accuracy 
up to the order 3PN. 

Our source modeling, as well as the one used in the EFT and ADM approaches, consists 
in representing each compact object as a (spinning) point particle whose internal structure 
is entirely parametrized by a set of effective multipole moments. The validity of this descrip¬ 
tion, which makes the calculations tractable analytically, relies on (i) the compact character 
of the bodies, and (ii) the weak influence of their internal dynamics to their “global” mo¬ 
tion in general relativity, often referred to as the effacement principle [49]. The foundations 
of this formalism were laid down in the seminal works of Mathisson [50 h 52| . Later Papa- 
petrou [53j, found a particularly simple form for the evolution equations (which comprise 
both the equations of motion and of spin precession) for dipolar particles, i.e. at linear order 
in spins. His derivation was improved and rephrased in the language of distribution theory 

f extensible beyond the dipolar model — 
55(. The dynamics of point particles with 


by Tulczyjew [54j], whose method — systematicall; 
has been recently applied at the quadrupolar level [ 
f ini te-size effects described by higher multipoles was thoroughly investigated by Dixon [56| 


59|, who constructed an appropriate stress-energy “skeleton” to encode information about 


the internal structure of the body while, on their side, Bailey & Israel proposed an elegant 
effective Lagrangian formulation (60|. Recently, Harte jgl] showed how the formalism of 
Dixon could be extended to self-gravitating systems, by constructing appropriate effective 
momenta and effective multipole moments evolving in some effective metric. 

In the present article, we are interested in the quadratic-in-spin contributions arising from 
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the quadrupolar moment of the compact object in the case where it is adiabatically induced 


by the spin |32|, [33|, [35|, [62|, as well as the simpler contributions coming from products 
of SO corrections. Because, in our source model, we replace extended bodies by point 
particles within a self-gravitating system, our approach must be regarded as an effective one 
and supplemented with some UV regularization procedure. A good choice is known to be 
dimensional regularization, with possible need of renormalization. We find however that, 
at this order, the so-called pure Hadamard-Schwartz prescription [63!] is sufficient, i.e. that 
dimensional regularization is not necessary. 

The paper is organized as follows. In Section [TT1 we explain how the dynamics of a 
test point particle endowed with a spin-induced quadrupolar structure moving in a curved 
background spacetime is described in the Dixon-Mathisson-Papapetrou formalism. We also 
write the equations of evolution for the particle worldline, as well as for the spin, under 
a convenient explicit form, and we define a spin vector of conserved Euclidian norm in 
terms of which our PN results shall be written. The validity of the model to describe 
the body dynamics in self-gravitating binaries is discussed. In Section IIIII dedicated to the 
computation of the next-to-leading order SS contributions to the PN equations of motion, we 
present expressions for the conserved energy in the center-of-mass frame, both for generic 
orbits and for the restricted case of circular orbits in the absence of precession. Finally, 
Section HVl sketches the derivation of the next-to-leading order SS contributions to the GW 
flux and includes a discussion of the impact of our newly derived terms on the phase evolution 
of non-precessing binaries in the frequency band of LIGO and Virgo. Because of the length 
of the equations, some results are relegated to appendices. Appendix [A] gives the explicit 
expressions for the relative acceleration and the precession vector in the center-of-mass frame, 
and Appendix [B] shows the relevant SS contributions to the source moments. We also give 
the explicit transformation between spin vector and spin tensor in Appendix 0 as well as 
the correspondence between our results and the ADM ones in Appendix 0 

We use the following conventions henceforth: 0{n) means 0{l/c n ), i.e. represents a 
contribution of the order (n/2)PN at least. Greek indices denote spacetime coordinates, 
i.e. /i = 0,1,2,3, while Latin indices are used for spatial coordinates, i.e. i = 1,2,3. Sym- 
metrization and anti-symmetrization are represented by, respectively, parenthesis and brack¬ 
ets around indices. We adopt the signature (—,+,+,+) and keep explicit both Newton’s 
constant G and the speed of light c. Finally the covariant derivative along the worldline is 
written as D/(cdr) = where is the four velocity of the particle, defined such that 


— 


vru 


= - 1 . 


II. DYNAMICS OF QUADRUPOLAR PARTICLES 

We shall now introduce the model we have adopted to represent the two spinning compact 
objects composing the binary as point particles. In Section III Al we display the Dixon- 
Mathisson-Papapetrou evolution equations for test bodies at quadrupolar order, set the 
covariant spin supplementary condition, and discuss its consequences. In Section III B1 we 
rewrite the equations of motion in terms of the 4-velocity and introduce a conserved mass. 
Section Hi Cl presents the construction of a spin vector with a conserved Euclidean norm and 
shows the precession equation it satisfies. Finally, Section Hi Dl explains to what extent the 
Dixon-Mathisson-Papapetrou dynamics can be used for the companions of a self-gravitating 
binary. 
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A. The Dixon-Mathisson-Papapetrou framework 


When describing the dynamics of a binary system of compact objects with masses 
A = 1,2, in the context of the post-Newtonian approximation, it is physically sound to 
model the two companions as point particles. Indeed, the ratio of the radii Ra ~ Gm^/c 2 
to the body separation r 12 is of the order GmA/(r 12 c 2 ), and thus much smaller than 1. The 
dynamics of test point-like objects including finite size effects has been investigated exten¬ 
sively by Dixon 156-59], who generalized the Mathisson-Papapetrou equations for spinning 
particles [.50, 51, HH, 64J by attaching arbitrary high-order moments to the individual bodies, 
beyond the monopole and the current dipole also referred to as the particle spin. It can also 
be derived from an effective Lagrangian-type approach for spinning particles, pioneered by 
Bailey & Israel 1I60|| (see also an extensive study for special relativity in (65j) and later im¬ 
plemented in EFT [21, 


, where higher-order moments appear as parametrizing couplings 
in the action to the value of the Riemann tensor and its derivatives on the worldline. 

The Dixon-Mathisson-Papapetrou equations of evolution for a spinning particle with 
quadrupolar structure read: 


D p a 
cdr 
D S’# 
c 2 dr 




(2.1a) 


(2.1b) 


where p a is the 4-momentum of the particle and u x = dx A /(cdr) the 4-velocity along the 
world-line. The anti-symmetric spin tensor S represents the effective 4-angular momentum 
of the object, while the (effective) mass and current type quadrupoles are encoded into the 
Dixon quadrupolar tensor Jp x p u ) which is only constrained at this stage to have the same 
symmetry properties as R pXpv . 

The stress-energy tensor T al3 of the model can be constructed after the Tulczyjew pro¬ 
cedure, by making the only assumption that its support is point like with at most two 
derivatives acting on the Dirac distributions, in the three following steps j55| : (i) write the 
most general symmetric tensor that involves up to two (covariant) derivatives of the particle 
scalar density 


n = 



cdr 


, 8*(x - y(r')) 

V~9 


( 2 . 2 ) 


where <5 4 (x — y(r)) is a 4-dimensional Dirac delta, with y[r) the particle worldline and x 
the field point; (ii) derive the hierarchy of equations verified by the coefficients of n in T pl/ 
due to the conservation equation \7 u T pi ' = 0; (iii) constrain those coefficients by solving 
all algebraic equations, which leaves two sets of ordinary differential equations. Identifying 
these two equations to Eqs. (12. ip yields the expression of T pi/ in terms of p^, S and J pA/w/ : 


T pv =n 


p^u^c + ^R (lJ Xrifr r )xpa c 2 


- V„ \nS p ^u y) ] - -V P V CT [nc 2 J p ^ u)ry ] 
o 


(2,3) 


It can be recovered with a smaller amount of calculation, further assuming that the system 
dynamics is governed by the effective Lagrangian of Bailey & Israel [60], by differentiating 
the resulting action with respect to the metric 
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As the spin tensor S p,v is anti-symmetric, it actually contains six degrees of freedom. 
Moreover, for an isolated body, the space-time components J° l of the total angular mo¬ 
mentum J = S pu /c in an appropriate asymptotically Minkowkian gauge represent the 
mass-type dipole of the object, and can thus always be taken to be zero. Similarly, for a 
test particle moving in a gravitational background, three degrees of freedom among those 
contained in the effective spin tensor are expected to be non-dynamical. They may be elim¬ 
inated by fixing the “center-of-body” reference point with the help of three independent 
space-time equations, globally referred to as the spin supplementary condition (SSC). The 
three remaining degrees of freedom correspond to the spatial components of the spin vector 
SA Various choices of SSC are possible (see for instanc e 166 ]i. Here we shall adopt, in 
keeping with previous works, the covariant (or Tulczyjew [54j) condition 

S^p u = 0 . (2.4) 


Assuming that the rotating bodies are always at equilibrium, we can reasonably expect 
their moments to depend on their masses, spins, as well as possible dimensionless parameters 
that characterize the internal structures. Notably, the spins may induce mass quadrupoles 
as they do for Kerr black holes. This effect produces spin square contributions that must 
be crucially taken into account at quadratic order in the spin variables. Tidal fields inside 
the bodies may also generate £ > 2 multipoles, but their leading order contribution to the 
acceleration would be ~ ( Ri, 2 /^ 12) 5 = 0 ( 10 ) for a compact binary, so that they can safely 
be neglected in the present work. 

As not all degrees of freedom in the Dixon quadrupole are physical, its value as a function 
of time cannot be uniquely determined by the internal dynamics of the body. In the adiabatic 
approximation, there exists a relation, valid along the particle worldline, between J px the 
4-velocity u p and the spin tensor S ,w . It can derived from an effective Lagrangian L$s built 
to be the most general Lagrangian — modulo perturbative redefinitions of the gravitational 
field, terms in the form of a total time derivative, terms that vanish under some given SSC, 
and 0(S 3 ) remainders — with the properties of: (i) being quadratic in S pi/ , (ii) depending 
on u p , the metric g jW , (derivatives of) the Riemann tensor, as well as some parameters 


characterizing the object [21, [33|, [35|. After redefining p p , S pi/ , we find that the stress- 


energy tensor associated with Lss coincides with that of Eq. (12.31) provided Jp x p v j s given 
by 


3k 

J pXpv = - S a[p u x] S} p u ^ ] , (2.5) 

rh c 4 

at any instant. The above expression properly describes the presence of a non-vanishing spin- 
induced quadrupole, with the source dependent constant k representing the quadrupolar 
polarisability. The mass parameter m is defined by p 2 = p a p a = —rh 2 c 2 . Notice that m is 
not a priori conserved. In fact, as shown below, its time derivative is quadratic in spin and 
cannot be consistently ignored at our accuracy level. 

The (contravariant) 4-momentum and 4-velocity of the particle are proportional when 
terms beyond linear order in the spins are neglected: p a = mcw“ + 0(S 2 ). Our first 
step will consist in expressing p a as a function of u a to quadratic order in the spin. We 
impose that the derivative along the worldine of the SSC (12.41) is zero, insert the equations 
of motion (12.lap and (I2.1bl) into the resulting identity, and use the fact that Jp x p v ~ 0(S 2 ) 
whereas S a/ 3 u f 3 ~ 0(S 3 ). This yields, at quadratic order in spin, 


p a = me u a 


Sap 

2 rh c 3 


-U T 


4c 


+ 0(S 3 


( 2 . 6 ) 
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We are now in position to write the spin evolution equation in a more explicit way. In 
the Lagrangian formalism, the effective linear and angular momenta are defined in a way 


that guarantees the conservation of the spin magnitude [42|, |62]. This conservation law is 


a remarkable feature of the spinning-particle dynamics. In our context, it will follow from 
Eq. (12.1|) for some class of supplementary conditions. In fact, it can indeed be derived ex¬ 
plicitly from those equations, for the form (12. 5 p of the quadrupole moment and the covariant 
SSC (12.4p . By substituting the 4-momentum (12.6p into equation (12.1bl) we get 


D5 , “' 3 

c 2 dr 


4c 

y l 


R'^J^ + u,uV>R\„r^ - u 1 uVr\^R\ iiv 


Oa]7 

- u^R^uV 1 —^- + 0(S 3 ) 


m c° 


(2.7) 


If we contract this expression with S a / 3 , we obtain S a pDS a/3 /(cdr) ~ 0(S 4: ) and, therefore, 
defining the spin magnitude as s 2 = S a pS a P/2 , 


ds 

dr 


(D(S 3 


( 2 . 8 ) 


This demonstrates that the spin magnitude is actually conserved at order 0(S 2 


B. Conserved mass and evolution equations 


Our next task is to investigate the issue of mass conservation at quadratic order in spins. 
For this purpose, let us compute the time derivative of the mass parameter rh. Using the 
equation of motion (12.lap and the Bianchi identities, we can write 

-fh^ = + 0 ( 5 =). ( 2 . 9 ) 

cdr cdr 6 

As the time dependence of Jp x p v is through the 4-velocity and the spin tensor, i.e. J pXtJ,u (r) = 
jp\jw ( m /3 ( r )^ ( r ))^ the fact that Did 3 /(cdr) ~ 0(S ) and DS' Q/3 /(cdr) ~ 0(S 2 ) implies 
the approximate conservation of the Dixon quadrupole: D J pXpu / (c&t) ~ 0(S 3 ). Now, 
substituting p with D/(cdr), we can write down the equation 


dr 


m - -R p x^J pXfl " 


0(S 3 ), 


( 2 . 10 ) 


which finally allows us to define a conserved quantity m as 

m = rh — ^R pXpv J pXpv . (2.11) 

Hereafter, the constant parameter m will be regarded as the effective mass of the particle. 
This mass is the one that appears in all our post-Newtonian results. By construction, it is 
conserved, like the spin magnitude. Substituting the expression (12.9p into Eq. (12.6p gives us 
the link between the 4-momentum and the 4-velocity: 

cap opu A 

P a = mcu a + -u a R pXpv J pXpv - + -J U ^ R X ^ J + 0 (S 3 ) > (2T2) 
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where m was just shown to be a constant parameter at order 0(S 2 ). We are then in position 
to rewrite the evolution equations for spinning particles to quadratic order in the spins, using 
the 4-velocity instead of the 4-momentum. Those are: 


D u a 
dr 


D S aP 

dr 


u‘ 0 DiW 1 pa >SZ_c, 

2 dr m 2 c 4 2 Xpu me 3 p x ‘ w m 
u a BR Xppv J xp ^ 4u B DR [a ppu J ,3]PIJ '' , , 0/c , 3 x 
o d r m 6 dr m 


4c 3 




- R lXpv u x u [ P 


^a]7 Ql 11 ' 


me 


+ e>(s 3 ). 




(2.13a) 


(2.13b) 


C. Definition of a spin vector and equation of precession 


From the anti-symmetric spin tensor S a ^, we define the spin 4-covector S a as 

S a = -\e aBpu —S lw , (2.14) 

2 me 

where e a / 3 pL , = \J—g rjap^ denotes the covariant Levi-Civita tensor, with g a / 3 pi , being the 
completely anti-symmetric symbol that verifies 770123 = A and where g = det g pv is the 
determinant of the metric tensor in generic coordinates. The tilde on this covariant spin 
vector will allow us to distinguish it from the Euclidean conserved-norm spin vector we shall 
introduce below. Notice that S a automatically satisfies S a p a = 0 and thus carries 3 degrees 
of freedom as required. If we contract the above equation with e" 7P °p 7 and use the SSC 
p u S pv = 0, we can invert Eq. 02.141) and obtain the spin tensor in terms of S a : 

S *P = PlL§ u + 0(S 3 ). (2.15) 

m c 

Remembering that pr = meu 13 + 0(S 2 ) at the linear-in-spin level, it is straightforward to 
check that S a S a = s 2 , by virtue of the relation e a/ 3 pu e xpar = —4 . 

To derive the evolution equation for the spin 4-covector, we differentiate Eq. 02.14j) with 
respect to the proper time, which yields 


Pg Q 

dr 




V n x ^- 

4 11 Act p U 


me 




+ 0(S 3 


(2.16) 


In what follows, we shall explicitly resort to our particular form (I2.5[) for ; relevant in 
the case of a spin-induced quadrupole. It will be convenient to investigate each term on the 
right-hand side of Eq. (12.161) individually. With our definition of S p , the first term there 
reads 


T^afipvR 


P 


s ap s 


Xap 


U 


2 Uk6 


crpK'Y Sry 


VA + o(s*). 


(2.17) 


4 ' ' ''“ r me 2 me 

At this stage, it is useful to introduce the gravitomagnetic part of the Bel decomposition of 
the Riemann tensor 




(2.18) 

























where *R is the self-dual Riemann tensor defined by 

*E> _ K 7 p 

g t a/3 ' 


(2.19) 


Physically, the tensor represents the tidal current-type quadrupole in the relativistic 
theory of tides. We can now put Eq. (12.16ft in the form 


1 P P ls xS ap S^ 1 e S y S p ^ 

^aP»vR Xap U mc ~ 2 H Ua mc +0 ( S >- ^ 2 - 20 ) 

Let us focus next on the second expression on the right-hand side of Eq. (12.16P . After 
substituting the value for J uXa P therein, we rewrite the resulting expression in terms of the 
gravitoelectric part of the Bel decomposition of the Riemann tensor 

= —RnXv P u x u p , ( 2 . 21 ) 

which is nothing but the tidal mass-type quadrupole generalizing that of Newtonian gravity 
(up to a factor 1/c 2 ). Next, we directly replace the spin tensors with their corresponding 
spin covectors in Eq. (12. 16ft . hence: 


^ - K e a ^G p ^ + 0(S 3 ). ( 2 . 22 ) 

dr 2 mc mc 

Finally, after setting 

Cl af3 = — [u [a H p] x - k e a ^ u u p G uX } , (2.23) 

771 C 

the spin precession equation for the covariant spin vector takes the form 

—- = n afi SP + 0(S 3 ). (2.24) 


The anti-symmetric tensor Q a p may be interpreted as a spin-precession frequency tensor. 



The passage to spin 3-vectors is achieved by introducing a direct orthonormal tetrad e£. 
The underlined index represents the vector label, which we may be viewed as the tetrad 
index, spacetime indices being represented by Greek letters and spatial indices by Latin 
letters as usual. Posing e 0 M = u p , we see that 

So = W = 0(S 3 ), (2.25) 


which means that Sq may be neglected. The squared Euclidean norm of Sa is then given by 

5^S b - = 7 = ~S,S p = s 2 , (2.26) 

with = g^ u + u^Uy. In words, the spin vector Sa has a conserved Euclidean norm. To 
define the spin variable uniquely in some coordinate grid, we still need to specify the choice 
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to be made for the spatial part of the tetrad. Considering that 5—e^e^,- = 7 ^, a natural 
choice is to take for the unique symmetric positive-definite square root (in the matrix 
sense) of 7 ^. The complete expression for the tetrad is 

(2.27) 

with v M == cu^/u° denoting the coordinate velocity. After projection on the basis vec¬ 
tors ( 12.27j) . the precession equation for the spin vector becomes 



d4 

dr 



where we have introduced the rotation coefficients for the tetrad 



De^ 

dr 


(2.28) 


(2.29) 


and where flaj 3 — e e ■ Now, as d/dr = u°d/dt, it is convenient to define an anti¬ 

symmetric precession frequency tensor associated with the coordinate time as 


^0 T --aBJ ■ 

Since S- is negligible, the precession equation reduces to 

dSi 


dt 


= n ij si+o(s- i ). 


(2.30) 


(2.31) 


Moreover, from the equality e 0 M = iC, it follows that the first term on the right-hand side 
of Eq. (12.231) vanishes when projected on spatial tetrad indices, so that 

n ij = -K£ ijk G^-^ : (2.32) 

- - me 

where £jj k or £ij k (indifferently) denote the Euclidean Levi-Civita symbol, with normalization 
£ 123 = £123 = 1, which is linked to the four dimensional Levi-Civita tensor by the relation 
t(H'ik ~eik‘ 

In the rest of the paper, we shall use a conserved Euclidean spin vector S with spatial 
components S l in harmonic coordinates such that 

S l = S-. (2.33) 


Because of the anti-symmetric character of fly, we can finally rewrite the precession equation 
in terms of a precession vector IT = —£^--Vtj k /2 as 


= e ijk WS k + Q(S 3 ). 


(2.34) 


It is the above precession vector IT that will be computed, along with the equation of motion, 
in Section m Our results will be displayed either in terms of the vector S l or of the spatial 
components of the spin tensor . 
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D. Application to self-gravitating binary systems 


Although the evolution equations (12.ip originally obtained by Dixon are only suitable to 
describe the dynamics of test particles, their rederivation based on the method of Tulczyjew 
or the Lagrangian approach of Bailey & Israel, regarded as effective field schemes, holds for 
self-gravitating N point-like body systems. Nonetheless, the validity of the point particle 
model breaks down at UV scales where the post-Newtonian expansion cannot be applied, 
i.e. for rA ~ Ra, with r being the distance between the particle representing the body A 
and the field point x. In particular, some infinities arise when computing the gravitational 
field iteratively due to divergences at the particle positions da■ The situation is even worse 
as we make x tend towards ija- 

As usual, those infinities are cured thanks to dimensional regularization, which preserves 
the invariance under diffeomorphism of general relativity, combined with some renormaliza¬ 
tion procedure. For an appropriate choice of the space dimension d, the field remains weak 
near = 0 and can be computed perturbatively in the post-Newtonian approximation. 
We are confident that this leads to the correct PN dynamics because: (i) the result for the 
acceleration is unambiguous up to the order 3.5PN for binaries of spinning compact objects, 
(ii) it is equivalent to that obtained from other methods (see the review paper (?;] for ref¬ 
erences), and notably from the approach a la Einstein-Infeld-Hoffmann used by Itoh (67) in 
the case of spinless bodies where no regularization is needed. Those cautions being taken, a 
self-gravitating system of N spinning bodies endowed with a quadrupolar structure may be 
modeled by means of the following effective stress-energy tensor, which generalizes that of 
Eq. (12.3j) : 


A=l,2 


n A IpMc+^V 


(m T v)\pa 2 


Xpa u A 


V p ( n A S&u) ) - ( n A c z J‘ 


„2 TpODo' 


(2.35) 


where the subscripts A indicate the particle label. 

The presence of poles oc e~ k in the metric at a given post-Newtonian order, with e = d — 3 
and k being a positive integer, may generate contributions in the source for the next or¬ 
der that could not be recovered by resorting to a purely three dimensional regularization. 
However, in the absence of such subtleties, the so-called pure Hadamard-Schwartz regular¬ 
ization [68j is sufficient to get the correct result. This prescription essentially relies on a 
specific use of the Hadamard partie finie regularization, which we shall briefly discuss now 


fiat)- 


(the reader will find more details in Ref. 

Let us consider a function F(x) with the same regularity properties as those arising in 
our problem, i.e smooth everywhere except at some singular points %/a (A = 1 , 2 ,..., A) in 
the neighborhood of which its admits an expansion of the form 


F ( x ) = r Afp( n A ) + °( r A) (2-36) 

po<p<P A 


for any integer P, with ua = (x — Ua)/pa■ Such a function is said to be of class T. Its 
Hadamard partie finie (F)a is then defined as the angular average of the finite part Afo( n A)'- 

(F)a = j , (2.37) 
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where df2^ denotes the elementary solid angle with direction tia centered on ija- The oper¬ 
ation of taking the Hadamard partie hnie is not distributive with respect to multiplication 
in the sense that, for another generic function G(x) of class J c , (F)a(G)a ^ ( FG)a in 
general. Moreover, it does not respect the Lorentz invariance. Because of the first of those 
two unpleasant features, the so-defined regularization is fundamentally ambiguous as such. 
Howbeit, it can still be used in practical computations provided it is supplemented by some 
additional prescription. In the PNISH approach, the post-Newtonian metric is constructed 
iteratively with the help of PN potentials. Those are elementary bricks satisfying a wave- 
type equation (more details are provided in Section IIIID . A convenient prescription is to 
define the value of a product FG of two potentials (or potential derivatives) evaluated at 
point yA as (F)a(G)a- Similarly, the regularized product of a potential F and an arbitrary 
smooth function a(x) will be given by a(yA)FA- 

Divergent integrals are cured by applying another kind of Hadamard partie hnie regular¬ 
ization. The regularized value of an integral with class-J 7 integrand is calculated in three 
main steps: (i) balls of radius rj centered on the singular points are extracted from the 
integration domain; (ii) terms that diverge near rj — 0 are removed; (iii) one goes to the 
limit rj —> 0. The singularities that generate poles in dimensional regularization produce 
logarithmic divergences in the Hadamard one. Those are associated with cutoff parameters 
Syi entering terms such as In(p/s^). For consistency between the two kinds of Hadamard 


regularizations, all derivatives must be evaluated in the sense of distributions 69]. The 


action of the three dimensional Dirac delta 5a = 5 3 (x — ua) on test functions must also be 
generalized to J r -class functions by posing F5a = ( F)a5a■ 

In this context, the pure Hadamard-Schwartz regularization is an ensemble of prescrip¬ 
tions designed to yield results that are “as close as possible” to those obtained through 
dimensional regularization. Those prescriptions demand: (i) to evaluate monomials of the 
form a(x)Fi...F n , where a(x) is a smooth function and the F^s are (derivatives of) T- 
class potentials, as ct(yA)(Fi)A---(F n )A] (ii) to evaluate divergent integrals by means of the 
Hadamard partie regularization for integrals; (iii) to extend the definition of 5a as explained 
above; (iv) to compute all derivatives in the sense of Schwartzian distributions. 

The absence of logarithmic cut-offs in the SS piece of the metric up to the order 3PN 
suggests that dimensional regularization may safely be swapped for the pure Hadamard- 
Schwartz one at this accuracy level. The insensitivity of the calculations to the choice 
of regularization procedure has been checked explicitly by evaluating source terms of the 
type FG5a in the stress-energy tensor as ( F)a(G)a5a , thus violating the pure Hadamard- 
Schwartz prescription. The results have always turned out to be unaffected by such modifi¬ 
cations. 


III. NEXT-TO-LEADING ORDER CONTRIBUTIONS TO THE POST-NEWTO- 
NIAN EVOLUTION 


We now turn to the computation of the dynamics of a binary system in the post- 
Newtonian approximation, at next-to-leading order for the quadratic-in-spin effects, i.e. 
at order 1/c 6 (or 3PN) in the equations of motion and at order 1/c 5 in the equations of 
precession. We will recover the results for the dynamics obtained in the ADM [62l . I70H73 


and EFT jH], 22, 3f| H, 39, 74j approaches, and extend them towards the completion of the 
calculation of the gravitational waves energy flux. 

We start with some general definitions in Section IIII Al Next, we introduce a set of 
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potentials parametrizing the PN metric in Section UlI Bl and express the quantities of interest 
in terms of these potentials. In Section IIII Cl we present their computation, and finally in 
Section HlI Dl the results obtained for the dynamics as well as various tests of their correctness. 
The lengthier calculations are all performed by means of the algebraic computing software 
Mathematica® supplemented by the tensor calculus package xAct [75l ] . 


General definitions 


The two objects are represented as quadrupolar point particles as explained above. An 
important ingredient of the formalism is the treatment of the infinite self-field of the point 
particles, essentially represented by means of Dirac deltas, through the pure Hadamard- 
Schwartz regularization procedure discussed in Section Hi Dl The distributional contributions 
yielded by derivatives are handled by using the Gel’fand-Shilov formula 76]. We found that 


at this order in spin, we have to keep track of distributional contributions in the metric 
itself to obtain the correct result for the wave generation formalism, as will be detailed in 
Section II V Bl 

The general structure of the equations of motion and precession is as follows: 


A - A n 4- —A 1 ™ 4- —A 2PN 4- J_4 2 - 5PN 4- ]_A 3 ™ 
71 — ^NS + 2 + c 4 ^NS + 5 ^i NS 4- 


- a n 

,3 ^SO 


4<o N + ®^s + 4< N + 0(7), 


«= 4«Ss +4«ns n + 4^0 +4«™+o(6), 


(3.1a) 

(3.1b) 


where the spin order in Eq. (I3.1b|) indicates the contribution in f2 itself, rather than in 
s = n x S (notably the SO terms feature the constants and actually correspond to SS 
terms in S). The 2.5PN NS terms in the acceleration are the first manifestation of radiation 
reaction. 

We use the same notations as in previous works. Three-dimensional indices are repre¬ 
sented with Latin letters a, &,... or i, j ,..., and are risen or lowered with the Euclidean metric 
®; we do not distinguish between upper and lower indices. We sometimes use boldface for 
Euclidean vectors. The positions and velocities of the two bodies are denoted by y[, y 2 and 
v\, v 2 . Apart from the separation distance ri 2 = 1 2 / 12 1 = \yi — 3 / 2 1 which we have already de¬ 
fined, we shall need the separation direction n\ 2 = {y\ — y l 2 )l r i 2 - The symbol 1 44- 2 indicates 
the same expression as the one before it, with the label of the two particles exchanged. The 
results are expressed in terms of the spatial components S \ : , S 2 of the spin tensor S^ v , as 
well as the spin vectors S\, S 2 of conserved Euclidean norm as defined above, in Section Hi Cl 
The mixed components S° l of the spin tensors can always be eliminated with the help of the 
spin supplementary condition (I2.4j) . We allow repeated indices in scalars quantities enclosed 
by parenthesis, in the absence of a risk of confusion. 

In harmonic (or DeDonder) gauge, the gravitational field equations can be rewritten as 


ah + a"" [ft] = Lhv" 


(3.2) 


where the stress-energy pseudo-tensor r ^ includes both matter and field contributions 
through and A MI/ , the latter source term being at least quadratic in . The field equa¬ 
tions (I3.2p . when iterated order by order, yield a solution expressed formally in terms of a 
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hierarchy of potentials of increasing complexity and post-Newtonian order (see Refs. 771.178 
for the precise definition of this iteration in the near-zone). 

Since on the one hand we are working at order 1/c 6 in the equations of motion, and 
on the other hand the spin contributions always come at relative 1/c order at least, only 
the so-called 2PN metric and potentials (i.e. necessary for the 2PN non-spinning case) are 
required. In fact, we will see below that, among the potentials arising at the order 2PN, 
only X turns out to be needed. For completeness, we quote here the result of the iteration 
for the 2PN metric, which reads 


9 oo = -1 + ^V-^V 2 + -Ax + VjVi + ^)+ 0(8 ), 


6 J 


90i = --7^ - ~rR t + 0(7) , 

9ij &ij 


2 2 
i + 4^ + 4 r 2 


+ Wij + 0(6) . 


The potentials therein are defined a£] 


V = D 




Vi = n£( 


-47 tG cr] 
-47tG (ji 


X = D 


-i 

n 


4v iG Van + W^V + 2V i d t d i V + Vdf V 




ZdiVAV, 




-471 tG (Vai - Va) - 2d k Vd,V k - -d t VdiV 


Wij = 04 [-4t rG {(Tij - Sijcj kk ) - diVdjV] , 
where the a, ai , aij quantities are convenient matter source densities defined as 


1 


a=^(T 00 + T n ), 


^ = ~T 0i 
c 


a i:i = T ij , 


(3.3a) 

(3.3b) 

(3.3c) 

(3.4a) 

(3.4b) 

(3.4c) 

(3.4d) 

(3.4e) 

(3.5) 


while stands for the PN-expanded retarded d’Alembertian operator acting on a function 
f(x, t) as [77, [78 


(□7 zf)(x,t) = ~^-J2~^ FFb=0 [ d?x'(\ x '\/ro) B \x 


— a;'! 71-1 f 




(3,6) 


n> 0 


Here FPb=o denotes the so-called Finite Part regularization, and r 0 is an associated arbitrary 
length scale. This regularization is used and described in Section IIV Al for the wave gener¬ 
ation formalism, but in Eq. 


it cures the divergences of the near-zone post-Newtonian 
metric at infinity rather than the divergences of the multipolar far-zone expansion at the 
origin. At the order we are considering here, it does not matter for the equations of motion. 
In particular the final results are all independent of the scale rjl. 

1 Possible contributions to the metric of non-linear tail terms, which are not made of (products of) elemen- 


77, 


78|. 


tary potentials defined by means of the operator , do not arise below the order 4PN 
2 At the 3PN non-spinning order, the scale r 0 does appear in the final results for the dynamics, but it 
disappears when considering gauge-invariant expressions such as E(ui), the conserved energy as a function 
of the orbital frequency. 
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B. Matter source and equations of motion in terms of potentials 


In this section, we introduce convenient additional definitions for the matter source in 
the PN context. In the covariant expression (j2.35jl . the worldline integration contained in 
the particle densities n A [see Eq. (j2.2jl ] can be performed explicitly in a definite coordinate 
grid (t,x). This results in 


T afS = Y, U A S A + + V„V„ 

A= 1,2 - 


(37) 


where we have defined 


rra/3 _ 

u A — 


/3) , I d ( “ 

. I Cp A U A -h % A I > 


u° A ^=gV yA “' A 1 3^ ax ^ 


TjaPii _ 

U A ~ „ n 


(a 

zU A o A , 


u aV~9 

1 f Z 7 p{ap)u\ 

Ua ~^n\~r jA )■ 


(3.8a) 

(3.8b) 

(3.8c) 


Here u° = 1 / —g A p v fl v u / c 2 , = (c, v*) (so that = u°v^/c), and the label index A on 

metric-dependent quantities means that they are to be regularized according to Hadamard 
regularization at the location of the particle A. In terms of partial derivatives, we have 


£ \TX l> /> A + -j=d r .{Tfn A ) + -±=d„ w (T; 


A= 1,2 


1 




(3.9) 


with 


q~ a P TT a P orb tjP)p v I (p> rr^)A/ii3 

’ A ~ U A + Apu U A + W ^)A U a 


A f.L , 

-,(a 




(TTP , _ , 

A pX\ w A Apv A 

-a/3_ /— -(tjoiPp . r /i TT a(3v\ or (a TT P)pv\^ 


Kn . 


= v/75(ty w + ry A t^ 


OpW TTh 

Zi Au\ u A 


Tf ^ 




(3.10a) 

(3.10b) 

(3.10c) 


By using Eqs. (13.81) . and the definitions of a, cq, oy,- given in Eqs. (13.5ft . we arrive at the 
following expressions in terms of metric potentials 


a = mi§i 


1 + ?(?* - v ') + ?(?4 - - 2,v + ~2^ + - 2 v 2 ) 


1 /7 


mic J 


2(Sf<<9 b H) - 4(S?d b V a ) ) + ^ ( - 8(f2 a <) + ^ 


1 ■>-. 1. 
/d 

11 


.6 


16 


33 


- 10vi(v“V a ) - 4(KK) + 2 (v“v1W ab ) - 3vfW - 8Z + —vfV 

O 

- 4«y a )y + 2WV + ^vfv 2 - \v* - 4X - AA_{S\ a Sfd ia V) 

4 O ZTH^ 




d k { 5i 


noka„,a oka 

^ + 4(-4VX + 4K, 
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rriic 6 

k 1 


sf sf a a y - (sf sf )<9 fe i/) 


+ + (■-SWV + 

«1 . f- W 1 , 3,^ 1 


miy^g 


dkl{ $1 


Cjka Cjla 

2 c 4 


+ ^(sMIj”? + 5 1 ") - 5 sf“s“<t.‘ 


+ i(5f 6 Sf)«M + Sft-SpDfof - Sj*«?®i)) 


+1^2+0 

(Ti = mi 5, 


-.7 




Hup 


2m!C 3 


+ 4f-^i - 4«KK - 2 ^i + + - V" 2 ^ 


2 c 3 x/=^ 
_ *h 

2mic 4 i/^^ 
«i 


d t Ui^X + 


mic 4 ^ 

+ 1^2 + 0^ 


(?ij = rri^i 


d> 


d t d k (s 1 s?s’r 

1 


<9fc 


2 2 

ik Qka rt .a rt .i x 


5i 


Qz/c Qka„.a n ,i 
£l_ °1 ^1^ 1 

2c 2c 3 


fc« 


c I ■*■ okanlcL^.i ^ nia ola n ,k ^ nia oka„,l 

Oi[-S 1 S 1 v 1 - -S 1 S 1 v 1 - -S 1 S 1 v x 


v i v i + — 


+ 


<V-0 


d k 


/>i(^si k vi+±si k v{ 


+1^2+0 


(3.11a) 


(3.11b) 


(3.11c) 


Here we have dropped the indices on the metric potentials, as we recall that, according to 
the pure Hadamard-Schwartz regularization, we can indifferently consider any quantity in 
factor of a Dirac delta as regularized, according to the rule FS\ = (F)i<5i. Notice however 
that the 1/^J—g prefactors must not be evaluated at point %)\. i.e. they are still functions of 
the field point x. These matter sources display explicit factors with spins and other without 
spin, but we should always keep in mind that there are secondary spin contributions coming 
from the potentials themselves. 

Let us now turn to the expression of the equations of evolution in terms of the metric 
potentials (13.41) . If we pose P ’ M = p^/m, using d/dr = u°d/dt, the covariant equations of 
motion (12. lap may be put in the form 


H P T) 1 

if = ■ F “ S ^ ’ < 3 ' 12 > 

where P /, = p^/m can be read from Eq. (12.12(1 . For a lower spatial index /i = i, we decompose 
P, as (coming back to notations for the body 1) 

p _ pNS , pSS 

bi - J- H + l \i , 
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(3.13) 


Fu = KT + F™ + Fi 


1 i 


Here, the order in spin refers to the order in the spin tensor as it reads in the formulas (12.12ft 
and (12.lap , but we should recall that there are also spin contributions coming from the 
potentials themselves, as well as from the replacement of accelerations using the equations 
of motion. We see from Eq. (12.12ft that Pi has no SO part in this sense, and its NS part 
comes from Pi = u°gi U v u /c. The NS part of F* comes from the usual connexion term in the 
geodesic equation, the first term in Eqs. (13.12j) . As the NS and SO parts can already be 
found e.g. in Eqs. (2.12) of Ref. T9j] and in Eqs. (3.7) of Ref. 24], we only display here the 
SS pieces: 


P (s 2 ) _ 


mfc 6 


Kl 


p(S 2 ) _ 


«1 


2 m\c A 


- SfSi b d t d a V + 2 S\ j S[ a d ja V b + S?(S?vZd ja V - 2St j vld jb V) 

}{s?s?a a v)*-s»s*M a v 

+ (SfSf)W + (S?S“ b )v“d ia V + S* j S*v$d ja V 
+ 2 SfS^djiV + SZ j S“\-2d jb Vi + 2d Tl V b ) 

Kl (sfsfyi^dAaV 


(3.14a) 


S?Sfd p d b V + 


,2^6 


mfc u 

^(sfsfjdfdy - ^(srsfd lb v)d,v + (sfsf)d a vd ia v 
- #sta.vw + s?s? 

+ 4 d a VdjiV) + STSf(2d t d J1 V b + -v\d ]lb V + \vd jib V 


+ SV{2 S{ H<^V-\ S fFA^V 

+ SfS h l j (2v a 1 d kia V j - 2v a 1 d kjl V a ) 


(3.14b) 


In these formulas, the potentials and their derivatives are to be understood as regularized 
at the location of body 1. 

For the equation of precession of the conserved-norm spin, we decompose similarly the 
precession vector into a NS and a SO part, before replacement of the potentials. We obtain 

(3-15) 

where (Q^) NS is given by Eq. (2.19) in j2^] and where 


( si i)so = ^b s ^ + 


rri\ c 


mic 0 
1 


v] (s°v l ;d ho v) - -(syi)d t d t v + sy;d ah v t 


- St(vW)diaV + - (S?V*) v“d ia V - {S“v“) d ia V a + S^diaVh 


K1 


- (S?d t d a V)v\ - -(S?vid ba V)v{ - (StvfidtdiV 
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(3.16) 


+ S[ (2 (yld t d a V) + (d a Vd 0 V) + (v“vid ba V) + <9 t V) 

- 3(S“d a V)d i V - ^(S a ] v°')v f ;d UJ V + S°(2d f d 0 V t + 2dAV a 
+ 2 vidbaVi + \v\d ia V - 3Vd ia V - 4v\d ia V b + 2v b 1 d ib V a ') \ 


The contributions featuring K\ come directly from the second term in Eq. (I2.30p . while the 
other contributions come from the first term there. The time derivatives of the velocities 
that enter the definition of the tetrad are replaced by the expression of the acceleration in 


terms of potentials, which include SO terms (as given for instance in Section 3 of 24|). 


C. Spin contributions in the metric potentials 


We now investigate the spin contributions to the metric potentials introduced in Sec¬ 
tion IIII Al As we are effectively working at the next-to-leading order, calculating these 
contributions from the results already presented above will be rather straightforward and 
we will only need to resort to well-known techniques. 

By inspection of the matter sources (13.111) . one can see that the SO and SS contributions 
to the metric potentials start at the following PN ordersll 


E so = 0(3 ), 
^ S ° = 0 ( 1 ), 
W so = 0 ( 1 ), 
Rf° = 0(1 ), 
Wg° = 0(3 ), 


E ss = 0(4 ), 
^ SS = 0(4 ), 
W ss = 0 ( 2 ), 
Rf = 0(4 ), 
Wg s = 0(4). 


(3.17) 


From Eqs. (j3. 14[) and (13. 16j) . we see that it is sufficient to compute the new SS contributions 
of the potentials V at the order 3PN (the leading order contribution at the order 2PN being 
already known, see e.g. Ref. 


33J), Vi at the order 2PN, and X at the order 1PN. 


We turn now to the caculation of E ss at the 3PN order which actually corresponds to 
the relative 1PN order. Truncating Eq. (13. 6 p appropriately, we have (dropping the FPb = 0 
regularization, which plays no role for compact sources) 


V ss = 


f / 

d 3 x' 

[J 

\x — x'\ 


-cr(f,x') 


1 d_ 

c di 


d 3 x'cr(t, x') + 


1 d 2 
2c 2 d t 2 


d 3 x'\x — x'\a(t, x') 


SS 


0 ( 8 ), 

(3.18) 

Because we are working at the next-to-leading order, various indirect contributions appear. 
Aside from the SS terms generated directly by the SS terms of a given in Eq. (13.1 lap , there 
are contributions from the SO 0.5PN part of V) in the SO 2.5PN part of a, from the SS 2PN 
part of V in the NS 2PN part of cr, and from the acceleration replacement featuring the SS 
2PN part of a 1 in the second time derivative of the integral of the NS Newtonian part of cr 
(in the third term above). 


3 It is implicitly understood that the orders n showed in the terms 0(n) below take their highest possible 
values. 
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In the following, to present the result in a more compact form, we adopt short cut 
notations: for any vectors a, b and spin tensors Sa,b, we define the scalars (ab) = a l b\ 
(, SaSb ) = S'aS 1 ^, ( aSAb ) = ahS^V, and ( ciSaSb b) = FS^S^b 3 (beware of the convention 
for the order of the indices on the spin tensors). 

We replace in Eq. (13.18[) the full expression (13.1 lap for er, perform integration by parts 
when derivatives of Dirac deltas occur, and compute the resulting integrals using Hadamard 
regularization, i.e. f d 3 xF6i = (F) l . The metric potentials can be considered as regularized 
when appearing in factor of a Dirac delta in the integrals, according to the pure Hadamard- 
Schwartz rule 63|] F5\ = (Fi)5i. 

An important point is that the derivatives have to be treated in a distributional sense. 
For the first time in our formalism, we have to take into account an essential distributional 
term in the potential V itself. The leading order result is indeed 


l/ ss = 


JA sfsfa ,Yl )+ i„ 2 + 0(6), 


(3,19) 


which, along with a non-distributional contribution, yields a distributional term given by 

(3.20) 


l/ ss — 

Mistr 




This distributional term will play no role in the derivation of the equations of motion them¬ 
selves, but it will produce a net contribution when computing the mass source quadrupole 
moment, as explained below in Section HV Bl Because, in this computation of the quadrupole 
moment, the V potential is only needed at the 2PN order, we will not need to consider pos¬ 
sible distributional terms at the higher 3PN order. 

Gathering the different non-distributional contributions we obtain 


K 


ss 

non distr 


GK\ 


2c 4 mir( 


(3(mSiSmO - (S'iS'i)) 


Gki 

C 6 771 i 


/ 3(niUi) 2 _ (vivi) _ 3Gm 2 {n 12 ni) 3Gm 2 {n l2 n 2 ) 3 Gm 2 

1 1 \ 4 r\ r\ 4rf 2 4 rf 2 2r 1 rf 2 


G 


3Gm 2 (ni 2 ni) Gm 


4 r\r\ 


12 


+ Gm 2 \ _ 3(u 1 S’iUi) s 
2 rfr 12 2r 2 rf 2 J 2 r\ 


+ (»,2 SAn„) 


15Gm 2 (ni 2 n 2 ) 9Gm 2 


4rf 2 


2 r 1 rf 2 


3 Gm 2 \ 
2 r 2 r\ 2 ) 



15(niUi) 2 3(uiUi) 3Gm 2 {n\ 2 ni) 3 Gm 2 

3 5 I V+2 1" 


4 r 1 


2rf 


2 rfr 


12 


2r\r\ 2i 


(v 1 S 1 S 1 v 1 ) 3Gm 2 {n l2 S 1 S 1 n 2 ) 

I o I 


+ W,S, n a ) ( 18G(ni2 ’ il) 


2r 4 

12 

6 G 


2r 4 

12 


.1 


3G(ni 2 SiS 2 n±) 


Zrf 2 


r i rf 2 J 

+ 102 + 0 ( 8 ). 

(3.21) 
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For the calculation of V) ss at its leading order 2PN, we proceed similarly as for V, but 
keeping only the first term in the expansion (j3.18[) . The calculation is simpler, with no 
indirect SS contributions. We get 

Vf S = (3(r.,S,S ini ) - (SiSO) v\ - + 1 «■ 2 + 0(6), (3.22) 

where we have included for completeness a distributional term completely analogous to the 
one discussed above for the potential G ss , but that will not contribute in the rest of our 
calculations. 

The computation of A ss is different, as it involves non-compact support terms. From 
Eqs. (13.41) , we see that the only 1PN SS contribution in X is 

X ss = D^ 1 [-2 diV^djV? 0 ] + 0(4), (3.23) 

with the leading order SO part of V) given by 

y. so = ^| + l«2 + 0(3). (3.24) 

We are working at leading order here, so that we keep only the first term in the expanded 
inverse d’Alembertian operator, which is just an inverse Laplacian. For the cross term, with 
derivatives of both 1 jr\ and 1 /r 2 , we use the function g = In [r\ + r 2 + r 12 ] which satisfies 
A g = l/(rir 2 ) (including the distributional part of the derivatives). With the notations 
d] = d/dy \, d 2 = djdy \, we can writcE 


A ’ 1 



dlA 



1 


rir 2 



(3.25) 


For the “self” terms, we can “factorize” the derivatives as explained in [81,] . Since we may 
ignore contributions of the form A -1 for £ + p even, we disregard possible distri¬ 

butional terms generated by space or time differentiation. After factorizing the derivatives, 
we transform them into derivatives with respect to y\ 2 and apply A -1 straightforwardly on 
the argument. The relevant formula is 
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d “ ( 9 du K 
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T (Sijfikl T T &il&jk) A j ^ 2 

A / J 


(3.26) 


Gathering these contributions, we find the following simple expression for the leading-order 
SS part of the potential X: 


X ss = 


G 2 


2c 2 rf [4 


-(SkS,) - (n 1 S 1 S 1 n 1 ) 


+1 ++ 2 + o(4), 


(3.27) 


where we keep the derivatives in the second term unexpanded. 


Including properly the regularization FP b —o would yield an additional constant contribution 


which 


would vanish after applying the derivatives. 
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D. Results for the evolution equations 


Using the results of the previous section and the NS and SO parts of the metric potentials 
that are already known, we are in position to complete the calculation of the equations of 
motion and precession (13. 14j) and (13 .161) . The results for the accelerations a 12 and the 
precession vectors 2 must pass several tests checking their validity. 

The first one is to make sure of the existence of a set of conserved quantities, in the 
absence of reaction reaction at this order, associated with the Poincare invariance of the 
problem: a conserved energy E, an angular momentum J, a linear momentum P, and a 
center-of-mass integral G. We were actually able to construct all those quantities explicitly 
by guess work. The higher-order terms in the precession equations intervene only in the 
conservation of the angular momentum, whereas the higher-order terms in the equations of 
motion intervene in all other conservation relations. We shall exhibit below the expression 
of the conserved energy, which will be later used to control the phase evolution of the binary 
in the case of circular orbits through the balance equation as explained in Section 11V Cl 
Another test consists in checking the Lorentz invariance of the dynamics, which must be 
manifest since the harmonic gauge choice is Lorentz-preserving. We use the same method 
as in Refs. 23, l24j . to which we refer the reader for more details, and find that our results 


pass this second test 

As the 3PN SS dynamics has been already investigated in both the EFT 


21, 22, 35, 38 


39], |74] and the ADM j(k| 7tJ [Z3[ approaches, we must be able to recover their results in our 
scheme. The equivalence between the ADM and EFT description has been shown to hold 
in Refs. jUl, E3, H] , so that we will only compare our results to the ADM ones, in keeping 
with our previous works. We present this comparison, and the resulting transformation 
from harmonic to ADM variables, in Appendix [D] The agreement with the ADM results 
also validates the test-mass limit of ours. 

Because the expressions produced are rather lengthy, we will give directly their reduced 
version in the center-of-mass (CM) frame. As in our previous works, this frame is defined 
as the one where the center-of-mass integral Gi (which is such that c\G/dt = P and hence 
d 2 G/dt 2 = 0) vanishes. We define x = rn = yi — y 2 the separation vector of the binary, 
v = dx/dt the relative velocity, m = m± +m 2 the total mass, v = mi r m 2 /m 2 the symmetric 
mass ration and 5 = (mi — m 2 )I'm the mass difference. We also use, for convenience, the 


same spin variables as in the previous works 23], |45j, namely 


S = Si + S 2 , S = m (— - . (3.28) 

\m 2 mi J 

The vectors S\ and S 2 are the conserved-norm vectors constructed in Section 111 Cl Addi¬ 
tionally, we will use the notation k + = K\ + k 2 and = K\ — k 2 . 

The positions of the two bodies the new frame read 

m 2 1 mi 1 . 

yi = —x + —z, y 2 = - -x + —z, (3.29) 

m cr m cr 

with 2 being a vector related to the center-of-mass integral G. In general, when working at 
the nPN order, only the (n — 1)PN expression of G (or z) is required. This can be checked 
explicitly from the Newtonian expressions in the general frame of the quantities of interest, 
as explained for instance in Ref. (25[. Thus, we would only need the SS 2PN expression of 
G in principle, but it turns out that there is no such contribution in G. We can therefore 
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translate our results to the CM frame using simply the same rules as in previous works: 
namely, we need the NS 1PN and the SO 1.5PN terms in z, as given in the Section 3 of 


Ref. 25 


For the SS contributions to the conserved energy, we find 


Ess = 


Gv 

c 4 r 3 


1 o 1 
7 e 4 S- 2 

4 <r 


Afi + 


1 Gm 


4 


(3.30) 


with 


e\ = S 2 (-2k+ - 4) + (SE) (-2 5k + - 45 + 2k_) + E 2 ((5k_ - «+) + v (2 k+ + 4)) 

+ ( nS ) 2 (6k + + 12) + (nE) 2 ((3 k+ — 35k_) + z/ (—6 k + — 12)) 

+ (n5)(nS) (65 k + + 125 — 6k_) , 
eg = S' 2 [(ra;) 2 ((65 k_ — 6 k + + 24) + v (6k + + 12)) 

+v 2 ((-2 5k_ + 8k + - 28) + i/ (2fi: + + 4))] 

+ (SE) [(nv) 2 ((—125k + + 485 + 12k_) + v (6 5k + + 125 — 30k_)) 

+v 2 ((105k + — 525 — 10re_) + v (25k + + 45 + 6k_))] 

+ E 2 [(nv) 2 ((65k_ - 6k+ + 24) + v (-95k_ + 21 k + - 72) + v 2 (-6k+ - 12)) 

+n 2 ((—55 k_ + 5k+ — 24) + v (5k_ — 11« + + 76) + z/ 2 (— 2k + — 4))] 

+ ( nS ) 2 [(m>) 2 zz (—30 k+ — 60) + v 2 ((60 — 18k+) + v (—6k + — 12))] 

+ (nS)(vS)(nv) ((—185k_ + 18k+ — 84) + v (12k+ + 24)) 

+ ( vS ) 2 (65k_ — 6k + + 28) 

+ (nE) 2 [(nv) 2 (y (155k_ — 15/c+) + v 2 (30k + + 60)) 

+u 2 ((95 k_ — 9k + + 48) + v (3 5k- + 15k + — 156) + u 2 (6k + + 12))] 

+ (nE)(vE)(nv) ((-185k_ + 18k+ - 72) + v (12 5k. - 48k+ + 228) + u 2 (-12 k+ - 24)) 
+ (uE) 2 ((65 k_ — 6k + + 24) + v (—65k_ + 18k + — 76)) 

+ (nS)(nE) [( nv) 2 u (—305k+ — 605 + 30/t_) 

+v 2 ((—185k + + 1085 + 18k_) + v (—6 5n + — 125 + 6k_))] 

+ (nT,)(vS)(nv) ((185k+ — 725 — 18k_) + v (65k + + 125 + 30k_)) 

+ (nS)(vT,)(nv) ((185/c+ — 845 — 18 k_) + v (65k + + 125 + 30k_)) 

+ (vS)(vE) ((—125k + + 525 + 12re_) + v (—24k_)) , 
eg = S 2 (—35k_ + 5k + + 8) + (SE) ((85k + + 85 — 8k_) + v (12re_)) 

+ E 2 ((4k+ — 45k_) + v (3 8k- — 11 k + — 10)) + ( nS ) 2 (9 5k- — 15k+ — 36) 

+ (nE) 2 ((12 8k- — 12k+) + v (—9 8k- + 33k + + 30)) 

+ (nS)(nE) ((—245 k + — 325 + 24k_) + v (—36k_)) . (3.31) 

The corresponding expressions for the relative acceleration a = a± — a -2 and the precession 
vectors S2 12 are provided in Appendix 1X1 

Finally, we further specialize our results to the case of circular, non-precessing orbits. As 
discussed in Ref. j33j, we have in fact three classes of orbits for the conservative dynamics. 
The CM expression are valid for general orbits, for which we make no assumption on the 
presence of precession and/or eccentricity. Quasi-circular precessing orbits correspond to 
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the case where we allow a generic orientation of the spins, but assume that the separation 
is constant at the SO level; as soon as SS and higher-order-in-spin terms are included the 


radius and orbital frequency become also variable on an orbital timescale. In Ref. 33] the 


definition of such orbits was investigated by perturbing orbital averaged quantities. The 
third and simplest class of orbits is that of the circular orbits with spins aligned with the 
orbital angular momentum and where precession is absent. As working at the next-to- 
leading order makes their determination more complicated, we leave the investigation of 
quasi-circular orbits for future work and focus on the circular, spin-aligned, non-precessing 
case. 

To present results for circular orbits, we use the same definitions as in previous works. 
We introduce a moving basis (n, A, £), with n denoting the unit vector along the separation 
vector, x = rn , £ = n x v/\n x v\ the normal to the orbital plane, and A completing the 
triad. When neglecting both radiation reaction and spin precession and assuming the spins 
aligned with £ the expressions for the relative velocity and acceleration become v = roj\ and 
a = — ruj 2 n, with u the orbital frequency defined by n = cnA. For the projected value of the 
(aligned or anti-aligned) spins along £ , we use the notation Sg = S ■ £. We also introduce 
the usual PN parameters 7 = Gm/rc 2 and x = (Gmoj /c 3 ) 2//3 , both of order 1PN. In the 
following, we only display the SS terms, and refer_the reader to Sections 9.3 and 11.3 of 
Ref. [?J for NS and SO contributions, and to Ref. j42j] for the newly computed cubic-in-spin 
contributions. 

First, we relate r to x by means of the equations of motion. We obtain the following SS 
terms for the PN generalization of Kepler’s law: 
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The result for the energy for circular, spin-aligned orbits is then 


+ 0(8) \ . (3.32) 
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(3.33) 


This expression can be shown to be in agreement, in the test-mass limit, with the energy 


of a test particle in circular equatorial orbits around a Kerr black hole 
control the phase evolution through the balance equation (see Section HVCD . 


82 . It is crucial to 


IV. NEXT-TO-LEADING ORDER CONTRIBUTIONS TO THE POST-NEWTO- 
NIAN GRAVITATIONAL WAVES ENERGY FLUX 

We now move to the computation of the 3PN spin-spin contribution to the energy flux 
radiated by the system. We start by briefly reviewing in Section IIV Al the basic elements 
of the wave generation formalism that we need here, before providing in Section ttV Bl some 
intermediate results useful in the calculation of the source multipole moments that are 
required to this order. The explicit results for the moments in the CM frame are relegated 
to Appendix [B] because of their length. Our explicit result for the GW flux is presented 
in Section IIV Cl for general orbits in the center of mass in the system and then reduced to 
the case of circular orbits in the configuration where the spins are aligned with the orbital 
angular momentum. 


A. Formalism 


We perform our calculation in the framework of the multipolar post-Newtonian approach 
to gravitational radiation. This formalism has been developed over many years, see e.g. 
[ 77 ; lEi • Since we will only use a simplified version of the full formalism, as we are 
working at next-to-leading order, we will refer the reader to [7] for a review, and give only 
a brief overview. 

The asymptotic waveform is defined from the transverse-tracefree (TT) projection of the 
metric perturbation, in a suitable radiative coordinate system X 11 = (cT, X), as its leading- 
order term in the 1/R expansion when the distance R = |X| to the source tends to infinity 
(keeping the retarded time Tr = T — R/c fixed). It can be parametrized using two sets of 
symmetric and trace-free (STF) radiative multipole moments, Ul of mass type and Vl of 
current type as 


h TT = 


4 G V£(n)£ 


c 2 R 


1=2 


U£\ 


UhL-2{Tr ) 


2 £ 


c(£+ 1) 


N m £mn(k Vi) n L— 2(-^~!r) 


ol w 


(4.1) 


where we denote by L — i\...if a multi-index composed of £ multipolar spatial indices A, ..., 
i^ ranging from 1 to 3. Similarly L — 1 = i\...if_\ and kL — 2 = ki\...if_ 2 ] Nl = N^.-.N^ 
is the product of £ spatial vectors N t . The transverse-traceless (TT) projection operator is 
denoted Vj^ = VikVji — -V;,'Pki where Vij = Sij — NiNj is the projector orthogonal to the 
unit direction N = X/i? of the radiative coordinate system. Like in the rest of this paper, the 
quantity is the Levi-Civita anti-symmetric symbol such that £123 = 1. The symmetric- 
trace-free (STF) projection is indicated using brackets or a hat. Thus Ul = Ul = U^l) and 
Vl = Vl = V(l) for STF moments. We denote time derivatives with a superscript (n). 
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In terms of these radiative moments, the energy flux into gravitational waves then reads 


+OO 


C=2 


G 


(l + l)(l + 2) ttWttW ,_ ^ + 2 ) _ v (1) V (1) 

(£-l)££\(2£+l)\\ L L c 2 (£-l)(£+l)\(2£+l)\\ L L 


(4.2) 


The Ul and Vj J can be expressed as (non-linear) functions of two sets of intermediate source 
rooted so-called canonical moments M L and Sl which are themselves related by a gauge 
transformation to a set of two so-called source multipole moments II, Jl (plus 4 gauge STF 
moments) which parametrize the most general solution to the Einstein equations outside 
the source. The differences between Ml and II (and similarly between Jl and Sl) arise at 
the 2.5PN order (see for instance [§}) and, since we are interested in SS effects which always 
add at least a factor 1/c 2 , we can safely ignore their differences. Using the same argument, 
we only need to consider the terms in the relation between the radiative moments and the 
canonical ones up to the order 2PN. Furthermore, we can also neglect the tail terms, which 
will only generate SS contributions at the order 3.5PN, so we finally have the simple relation 

(tf(j)ss = (4 2> )ss +0(7), (4.3) 

(Uj)ss = (4 ?> )ss + O (7) . (4,4) 

(U m ) ss = (42) ss + 0(7). (4.5) 

Noticing additionally that the leading order spin-spin contribution to any of the II or Jl 
(and their time derivatives) is of the order 2PN (as will be clear from the expressions in the 
next section), we can express the spin-spin flux in terms of the relevant source moments as 


•Fss 



1 A 3) t(3) 
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1 r(4) r(4) 16 (3) j( 3) 

189 45 v 



J_ j( 4) t(4) 
§4 J ijk J ijk 


+ 0(7) , 
ss 


(4.6) 


which requires computing the SS parts of I t j to the order 3PN and of J t j and Iijk to the 
order 2PN. We also need the NS parts of I t j up to the order 1PN and of J t] and I Vj k at 
the Newtonian order, as well as the SO contributions in l l3 and Jij up to the order 1.5PN 
and the leading 0.5PN SO contribution to J^, all of which are known from previous works. 
Remember that the spin-orbit contributions to mass (resp. current) type moments start 
at 1.5PN (resp. 0.5PN) order, and that time derivatives of non-spinning (resp. spin-orbit) 
expressions generate spin-spin contributions with an additional order 2PN (resp. 1.5PN) at 
least. 

The matching procedure at the core of the formalism finally allows us to express the 
source moments as closed-form integrals over space j87j. Instead of reproducing here the 
general expressions which can be found in Eq. (123) of Ref. [?J, we directly display below 
the terms that contribute to the spin-spin corrections at the required orders. They read 

('«>« = E /d 3 * 

(Jij) ss = |P J d 3 x e ah<u x l>( ff + O ^ , (4.7b) 

(Ii jk )ss = FP J d 3 x (-0 % fc E SS + O (i) , (4.7c) 


E + 
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14c 2 
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21c 2 qij q 


ss 
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(4.7a) 
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where FP# = o denotes a finite part operation defined by analytic continuation in the complex 
plane for the parameter B, which deals here with the infrared divergences at infinity. An 
arbitrary scale tq is introduced, which will play no role in the present calculation and has to 
disappear from gauge-invariant results. The basic “building blocks” £, £, and + t] entering 
the integrands are defined as 

_00 I 

£ = —- , £* = — , S ij = r ij , (4.8) 

where r ^ has been defined in (13.2p . and the overline indicates a post-Newtonian (near-zone) 
expansion. In identifying the relevant terms in (I4.7a| . we slightly anticipated on the results 
of the next subsection (see Eq. (I4.9a| ) and used the fact that the SS contributions to E, £,, 
and Ejj all start at the 2PN order at least. 


B. Computation of the source moments 

To obtain the relevant SS contributions to the source moments, we first express the sources 
E, E i and E^ in terms of the potentials parametrizing the metric and the matter sources a, 
ai and aij defined in (13.51) (the complete relations can be found, generalized to d, dimensions, 
in @1)- Taking into account the order of the spin corrections in these quantities, the only 
terms that yield spin-spin contributions to the orders we are interested in are 
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7 tGc 2 


SS 


(4.9a) 

(4.9b) 

(4.9c) 


The i nteg rals in Eq. (I4.7al) can now be performed using the standard techniques described 
in 6$, 691, handling the UV divergences of the integral through the Hadamard regularization 
and the IR divergences through the finite part operation FP^ =0 . 

We highlight here that the distributional parts of the sources have to be treated with 
care. In particular, for the first time, we encountered the situation where such contributions 
in the metric itself (more precisely in the potential V), and not just those coming from 
derivatives applied to the metric, have to be crucially taken into account. 

More specifically, the spin-spin leading order contribution in the potential V was com¬ 
puted in Eqs. (13.19j) and (13.20)1 and contains a term proportional to <5i which has to be 
accounted for when integrating the diVdiV term of (I4.9a.p in (14.7a In order to illustrate 
this further, let us focus on the second and third terms in Ess 


=ss 4V 
E v 


7 tGc 2 


diVdiV, 


(4.10) 


which we can rewrite using the identity ‘IdiAdiB = A (AB) — AAB — BAA, and the fact 
that AV = —47 tGct at leading order, as 


^ss 


1 1 


^ ■ 


(4.11) 
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By reinjecting this second form into f!4.7aH . integrating by parts, using Ax tl = 0 and treating 

the surface terms as explained in the Section IV D of @1, we readily see that E^ S actually 
gives a vanishing contribution to iy. If on the other hand one uses (I4.10p without including 
the distributional part of V, one obtains an incorrect non-zero result. 

Our explicit results for the SS contributions to the source moments reduced to the center 
of mass are presented in Appendix [B] 


C. Gravitational waves energy flux 

Using equation (14. 6 1) . our results for the source moments and the equations of motion and 
precession obtained in Section IIIIDI to compute time derivatives, we can finally compute 
explicitly the gravitational wave flux. We will give the result already reduced in the center- 
of-mass frame, and we use the same notations as already introduced in Section IIII Dl We 
obtain 

G 3 m 2 n 2 
‘ Fss = 5c 9 r 6 

with 

fl = S 2 [(nu) 2 (-312 ac+ - 624) + v 2 (288k + + 576)] 

+ (SE) [(nv ) 2 (-3125k + - 6245 + 312k_) + v 2 (2885k+ + 5765 - 288k_)] 

+ E 2 [(m>) 2 ((1565 k_ - 156k + + 18) + v (312k+ + 624)) 

+v 2 ((-1445/c- + 144 k+ + 6) + v (-288k+ - 576))] 

+ (nS) 2 [(nv) 2 (1632k+ + 3264) + v 2 (-1008k + - 2016)] 

+ (nS)(vS)(nv) (—696k + — 1392) + ( vS ) 2 (144k + + 288) 

+ (nE) 2 [(nv) 2 ((-8165k_ + 816k + + 18) + v (-1632k+ - 3264)) 

+v 2 ((5045k_ - 504k + ) + v (1008k+ + 2016))] 

+ (nE)(nE)(nn) ((3485 k_ — 348k + — 12) + v (696k + + 1392)) 

+ (uE) 2 ((—725 k_ + 72k + + 2) + v (-144k+ - 288)) 

+ (nS)(riE) [(nv) 2 (16325k+ + 32645 - 1632k_) + v 2 (-10085k+ - 20165 + 1008k_)] 

+ (nE)(vS)(nv) (—3485k + — 6965 + 348k_) + (nS)(vT,)(nv) (—3485k+ — 6965 + 348k_) 
+ (vS)(v E) (1445k + + 2885 - 144k_) , 

/ 6 ° = S 2 [(nv) 4 ((22745k_ + 12918k+ + 35436) + v (-14112k+ - 28224)) 

+(nv) 2 v 2 ((-25925k_ - 17544k+ - 51984) + v (17928k+ + 35856)) 

+v 4 ((3665k_ + 5034k + + 18276) + v (-4584k + - 9168))] 

+ (SS) [(nv) 4 ((106445k + + 506525 — 10644k_) 

+u (-141125k+ - 282245 + 5016k_)) 

+(nv) 2 v 2 ((—149525k + - 696725 + 14952k_) 

+v (179285k+ + 358565 - 7560k_)) 

+v 4 ((46685k + + 208125 - 4668k_) + v (-45845k+ - 91685 + 3120k_))] 

+ E 2 [(nv) 4 ((—53225k_ + 5322k+ + 9714) + v (47825k_ - 15426k+ - 64788) 


-f°-\ _ — 

3 U 21c 2 


o Gm i G 2 m 2 2 

„ ^6 "T __9 J 6 


(4.12) 
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+z/ 2 (14112/C+ + 28224)) 

+(nv) 2 v 2 ((74765k_ - 7476k+ - 14286) 

+/y (—6372<5 ac_ + 21324k + + 86316) + v 2 (-17928 k + - 35856)) 

+v 4 ((—23345 k. + 2334k+ + 3796) + v (1926 5k. - 6594k + - 23336) 

+z/ 2 (4584 k+ + 9168))] 

+ (nS) 2 [(nv) A ((129305k_ - 90570*;+ - 81570) + v (71520*;+ + 143040)) 

+(nv) 2 v 2 ((-81245k_ + 81636*;+ + 65220) + v (-62976*;+ - 125952)) 

W ((5705*;_ - 14778*;+ - 6546) + v (13632*;+ + 27264))] 

+ (nS)(vS) [( nv ) 3 ((-19752 5k. + 51816*;+ + 16464) + v (-29184*;+ - 58368)) 

+(nv)v 2 ((9522 5k. - 19890*;+ - 180) + v (4092*;+ + 8184))] 

+ (vS) 2 [(nv) 2 ((6378 5k. - 9114*;+ + 7794) + v (5100*;+ + 10200)) 

+v 2 ((-16685k_ - 324*;+ - 6478) + v (120*;+ + 240))] 

+ (n£) 2 [(?nvf ((517505*;_ - 51750k+ + 18420) 

+v (-486905k_ + 152190^+ + 60960) + u 2 (-71520*;+ - 143040)) 
+(nv) 2 v 2 ((-448805*;. + 44880*;+ - 8112) 

+u (396125*;. - 129372*;+ - 79608) + u 2 (62976*;+ + 125952)) 

+v 4 ((7674 5k_ - 7674*;+ + 3090) + v (—73865«:_ + 22734*;+ + 10884) 

+z/ 2 (-13632*;+- 27264))] 

+ (n£)(w£) [(nv) 3 ((—35784<5«:_ + 35784*;+ - 33534) 

+v (343445 k. - 105912*;+ + 48858) + v 2 (29184*;+ + 58368)) 

+(nv)v 2 ((147065k. - 14706*;+ + 7782) + v (-115685*;. + 40980*;+ - 7794) 
+z/ 2 (—4092*;+ — 8184))] 

+ (vT) 2 [(nv) 2 ((7746 5k. - 7746 k + + 14124) + u (-8928 5k. + 24420*;+ - 36432) 

+z/ 2 (—5100*;+ - 10200)) 

+v 2 ((-672 5k. + 672 k + - 2242) + v (16085k. - 2952k+ + 9788) 

+z/ 2 (—120 k+ — 240))] 

+ (nS)(riE) [(nv) 4 ((-1035005k+ - 709205 + 103500k.) 

+u (71520 5k + + 1430405 - 123240k.)) 

+(nv) 2 v 2 ((897605k+ + 718085 - 89760k.) 

+v (-62976 5k+ - 1259525 + 95472k.)) 

+v 4 ((—153485k+ - 86645 + 15348k.) + v (136325k+ + 272645 - 15912k.))] 

+ (riE)(vS) [(nv) 3 ((357845k+ — 154025 — 35784k.) 

+u (-14592 5k + - 291845 + 54096k.)) 

+(nv)v 2 ((—147065k+ + 81905 + 14706k.) + v (20465k+ + 40925 - 21090k.))] 
+ (nS)(vY,) [(nv) 3 ((357845k+ — 2405 — 35784k.) 

+u (-14592 5k+ - 291845 + 54096k.)) 


+(nv)v 2 ((—147065k+ - 51245 + 14706k_) + u (20465k+ + 40925 - 21090k_))] 
+ (vS')(vS) [( nv) 2 ((—154925k + + 230525 + 15492k_) 

+v (51005k+ + 102005 - 30612k_)) 

+v 2 ((13445k+ - 81885 - 1344k_) + v (1205k + + 2405 + 6552k_))] 
fl = S 2 [(nv) 2 ((—2772 5k_ + 23844^+ + 48872) + v (-1320k+ - 2640)) 

+v 2 ((25725k_ - 21028k+ - 41832) + v (720k+ + 1440))] 

+ (SE) [(nv) 2 ((266165k+ + 551765 - 26616k_) + v (-13205^+ - 26405 + 12408k_)) 
+v 2 ((-236005*;+ - 488725 + 23600k_) + v (7205*;+ + 14405 - 11008k_))] 

+ E 2 [(nv) 2 ((—133085*;_ + 13308*;+ + 1208) + v (34325*;_ - 30048*;+ - 61736) 

W (1320*;+ + 2640)) 

W ((118005*;_ - 11800*;+ - 4408) + v (-29325*;_ + 26532*;+ + 56288) 

+v 2 (—720*;+ — 1440))] 

+ (nS) 2 [(■ nv) 2 ((287885*;_ - 135588*;+ - 300528) + v (2028k+ + 4056)) 

+v 2 ((—123805*;_ + 77736k+ + 182752) + v (-1992*;+ - 3984))] 

+ (nS)(vS)(nv) ((-204725*;_ + 64056*;+ + 123000) + u (1932*;+ + 3864)) 

+ (vS) 2 ((46645*;_ - 14652*;+ - 27240) + v (-168*;+ - 336)) 

+ (nE) 2 [(nv) 2 ((821885*;. - 82188*;+ - 5604) + v (-298025k. + 194178*;+ + 264672) 
+v 2 (—2028*;+ - 4056)) 

+v 2 ((—450585k. + 45058k+ + 9700) + v (133765k. - 103492k+ - 185532) 
W (1992k+ + 3984))] 

+ (nE)(vE)(nv) ((-422645 k_ + 42264k+ - 9808) + v (195065k. - 104034k+ - 69804) 
W (-1932/C+ - 3864)) 

+ (vE) 2 ((96585 k. - 9658k+ + 4784) + v (-45805k. + 23896k+ + 9808) 

W (168k+ + 336)) 

+ (nS)(nE) [(nv) 2 ((-1643765k+ - 2824445 + 164376k.) 

+v (20285k+ + 40565 - 117180k.)) 

W ((901165k+ + 1837165 - 90116k.) + v (-19925k+ - 39845 + 51512k.))] 

+ (nE)(vS)(m;) ((422645 k+ + 464325 - 42264k.) + v (9665k+ + 19325 + 39978k.)) 

+ (nS)(vE)(nv) ((422645k+ + 507445 - 42264k.) + v (9665k+ + 19325 + 39978k.)) 

+ (v^)(wE) ((-193165 k+ - 184805 + 19316k.) + v (-1685k+ - 3365 - 18488k.)) , 
fl = S 2 ((165/c_ + 144 k+ + 368) + v (-576k+ - 1152)) 

+ (SE) ((1285k+ + 2245 - 128k_) + v (-5765k+ - 11525 + 512k_)) 

+ E 2 ((—645k_ + 64k+ + 24) + v (2725k_ - 400k+ - 384) + v 2 (576k+ + 1152)) 

+ (nS) 2 ((—485k_ - 432k+ - 936) + v (1728k+ + 3456)) 

+ (nE) 2 ((1925 k. - 192k+ - 16) + v (-8165k. + 1200k+ + 928) 

W (-1728k+ - 3456)) 
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+ (nS)(nL) ((—3845k + - 7845 + 384k_) + v (17285k + + 34565 - 1536k_)) . (4.13) 


After reduction to the case of spin-aligned, circular orbits, using the notations already 
introduced in Section 1111 Dl for the energy, we obtain 
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(4.14) 


Using this result as well as the expression of the orbital energy (13.33p . we can write the 
balance equation T = — dE/dt for circular orbits to obtain the phase evolution of the 
binary. Different ways of mixing analytical and numerical integration give rise to different 
approximants (see for instance [89j] for a comparison of these different approximants). For 
simplicity, we will give here only the phasing formula for the TaylorT2 approximant: we 
re-expand d(p = 2u>dt = 2u(—J 7 /(dE/dt)) and integrate term by term to obtain the phase 
of the wave (j) (here (j) is the phase of the leading 22 mode, hence the factor of 2) as a function 
of u> or equivalently of x. We get for the SS contributions 


(</>)ss = 
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The known NS and SO contributions are summarized in Sections 9.3 and 11.3 of [7j, and 


LIGO/Virgo 

10 M q + 1.4M 0 

IOMq + 10 Mq 

Newtonian 

3558.9 

598.8 

1PN 

212.4 

59.1 

1.5PN 

—180.9 + 114.0xi + 11.7y2 

-51.2 + 16.0x1 + 16.0x2 

2PN 

9.8 - 10.5 X ? - 2.9 X iX2 

4.0-l.lx?-2.2xiX2-1-1x1 

2.5PN 

-20.0 + 33.8x1 + 2.9x2 

-7.1+ 5.7xi + 5.7x2 

3PN 

2.3 - 13.2x1 - 1.3X2 
— 1-2X? - 0.2xiX2 

2.2 - 2.6xi - 2 . 6 x 2 
-0.lxf-0.2xiX2-0.lxl 

3.5PN 

-1.8+ 11 . 1 x 1 + 0 . 8 x 2 + (SS) 
—0.7xi - 0.3xiX2 

-0.8 + 1.7xi + 1.7x2 + (SS) 
-0.05xf - 0.2x|x2 - 0.2xixi - 0.05x1 

4PN 

(NS) - 8.0X1 - 0.7x2 + (SS) 

(NS) - 1.5xi - 1.5x2 + (SS) 


TABLE I. Number of cycles associated to the different PN terms in the phasing formula, between 
the starting frequency for advanced detectors (10Hz) and a cut-off chosen to be the Scwarzschild 
ISCO x = 1/6. We show the result for typical black hole/neutron star and black hole/black hole 
systems. Spin-aligned, circular orbits are assumed, and we use the dimensionless spins xa such that 
Sa£ = Gm? A XA■ We ignore contributions that are at least quadratic in the spin of the neutron star. 
We gather all contributions known to date, the ones still unknown are indicated in parenthesis. 


additional cubic-in-spin 3.5PN contributions can be found in 42]|. We give in table [J the 


number of cycles of the signal resulting from each term in the phasing formula, for the 
frequency band of advanced LIGO/Virgo detectors. Notice however that these results are 
illustrative, as they are specific to the TaylorT2 approximant and as these number of cycles 
give only a rough idea of the relevance of these terms in actual data analysis applications. 

We can check that our result (I4.14p is in agreement, in limit of a test particle orbiting a 
Kerr black hole, with the result of |90| obtained in the framework of black hole perturbation 
theory. We leave for future work the comparison of our results with the so far incomplete 
results (given only at the level of the multipole moments) of [43, |44], Natural extensions 


of the work presented here include the investigation of quasi-circular precessing orbits, the 
computation of the spherical harmonic decomposition of the waveform (or, equivalently, the 
full polarizations /i +iX ), and the implementation of these results for the factorized waveforms 


of the Effective-One-Body formalism with spins (see e.g. fill 194 
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Appendix A: Explicit results for the equations of evolution 


We gather in this appendix explicit results for the equations of motion and precession 
obtained in Section fill Dl that are too long to be shown in the main text. We present them 
already reduced to the center-of-mass frame. 

For the precession vectors Q) 2 , we have found simpler to keep the variables S\ t2 and k 12 
instead of S, £ and k±. We get 
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(A2) 


For the relative acceleration a 1 = — a 2 , we obtain (coming back to the S,H and 

variables) 


Mbs = 
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c 4 r 4 m 


1 i 1 

a d n + 

c 


4 " 4,0 1 ^ 2 


1 , 1 Gm 

F a 6,o + 


4 r “ 6a 


0 ( 8 ), 


(A3) 


with 


< 0 = S l 


(nS) (— 12k + — 24) + (n£) (—65 k + — 125 + 6k_) 

(nS 1 ) (—65 k + — 125 + 6k_) + (n£) ((65 k_ — 6k + ) + v (12k+ + 24)) 
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n" S 2 (— 6 k + 
(SE) (- 6 Sk + 


12 ) + E 2 ((35k_ - 3k+) + v ( 6 k+ + 12 )) + (nS) 2 (30k+ + 60) 
125 + 6 k_) + (nE ) 2 ((15k + — 155k_) + v (—30k + — 60)) 


+ (nS')(nS) (305 k + + 605 — 30k_) , 


a, 


6,0 


= S i 


(■ nS ) ((nv ) 2 ((60k + — 605k_) + v (60k + + 120)) 


+v 2 ((126k- - 36k + - 48) + v (-72k+ - 144))) 

+ (nE) ((ra;) 2 ((605n + — 1205 — 60n_) + v (305n + + 605 + 90n_)) 
+v 2 ((24k_ - 245n+) + v (—365n+ - 725 + 12/c_))) 

+ (vS)(nv) ((305n_ — 30n + + 84) + v (— 12 k + — 24)) 


+ (vE)(nv) ((—305n+ + 1325 + 30n_) + u (—65n + — 125 — 54n_)) 

+ E* ( nS)(nv ) 2 ((605n + — 60n_) + v (305n+ + 605 + 90n_)) 

+ (nE)(m ;) 2 ((—605n_ + 60n+ - 120) + v (305n_ - 150«+ + 240) + v 2 (-60n+ - 120)) 
+ (■ nv)(vS ) ((—305n + + 485 + 30n_) + v (— 6 6k + — 125 — 54n_)) 

+ (nu)(uE) ((305n_ - 30n+ + 96) + v (-245n_ + 84n+ - 276) + v 2 (12«+ + 24)) 

+ ( nS)v 2 ((—245n + — 245 + 24n_) + v (—365n + — 725 + 12re_)) 

+ (nE)v 2 ((245n_ - 24 k + + 24) + v (245n_ + 24 k + ) + u 2 (72 k+ + 144)) 

nS) 2 (nvfv (-210/C+ - 420) + (nv) 2 S 2 ((305n_ - 30n+ + 120) + v (30n+ + 60)) 

(nS)(nE)(nv) 2 v (—2105n + — 4205 + 210n_) + ( vS ) 2 (65n_ — 6 n + + 84) 

(nE) 2 (ra ;) 2 (u (1055n_ — 105n + ) + u 2 (210n + + 420)) 

(nv) 2 (SE) ((—605n + + 2405 + 60n_) + v (305n + + 605 — 150n_)) 


+ n 


+ (nv) 2 E 2 ((305n_ - 30n+ + 120) + v (-455n_ + 105n+ - 360) + v 2 (-30n+ - 60)) 
+ (nS)(nv)(vS) ((—305n_ + 30n + — 420) + v (60n + + 120)) 

+ (nE)(nv)(vS) ((305n + — 2405 — 30n_) + v (305n + + 605 + 30n_)) 

+ (nS)(nv)(v E) ((305n + — 4205 — 30n_) + v (305n + + 605 + 30n_)) 

+ (nE)(nv)(vE) ((-305n_ + 30n+ - 240) + v (900 - 60n+) + v 2 (-60n+ - 120)) 

+ (vS)(v E) ((—125n + + 1325 + 12 «_) + v (—24n_)) 

+ (v E ) 2 ((65n_ — 6 n + + 48) + v (—65«_ + 18n + — 180)) 

+ (nS) 2 v 2 ((60n+ + 120) + v (180n+ + 360)) 

+ S 2 v 2 ((— 6 6k- — 6 k + — 48) + v (—36n + — 72)) 

+ ( nS)(nE)v 2 ((605n + + 1205 — 60n_) + v (1805n + + 3605 — 180n_)) 

+ (nE) V ((30k + - 305n_) + v (-905n_ + 30n + - 120 ) + v 2 (-180k+ - 360)) 

+ (SE)u 2 ((-725) + i/ (—365n+ - 725 + 60n_)) 


+ E V (-24 + v (245n_ - 24n+ + 96) + v 2 (36n+ + 72)) 


+ v l 


( nS) 2 (nv ) ((—240k + ) + v (120k + + 240)) 
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+ ( nv)S 2 ((—12 5k,- + 60k + — 48) + v (— 24k + — 48)) 

+ (nS)(nE)(nv) ((-2405k+ + 2405 + 240«_) + v (120 8k+ + 2405 - 120«_)) 

+ (nE) 2 (m;) ((1205k_ - 120k+ + 240) + v (-605k_ + 300k+ - 480) 

+u 2 (—120 k + - 240)) 

+ (nv)(SE) ((72 5k + - 1445 - 72«_) + v (-245^+ - 485 + 72k_)) 

+ [nv)Y? ((—365k_ + 36k+ - 96) + v (245k_ - 96k+ + 240) + v 2 (24k+ + 48)) 
+ (nS)(vS) ((65k_ + 90k + + 84) + v (—36k+ — 72)) 

+ (nE)(vS) ((425k + — 42k_) + v (— 185k + — 365 + 6 /c_)) 

+ (nS)(vE) ((425k + + 365 — 42k_) + u (— 185k + — 365 + 6k_)) 


+ (nE)(vE) ((—425 k_ + 42k+ - 48) + v (125k_ - 96k + + 12 ) + v 2 (36k+ + 72)) 
a\ v = S i |"(nS') ((— 245k_ + 72k+ + 164) + v (36k+ + 72)) 

+ (nE) ((485 k + + 725 — 48k_) + v (185k + + 365 + 30k_)) 

[nS) ((485k + + 845 — 48k_) + v (185k + + 365 + 30k_)) 

+ (nE) ((48 k + — 485k_) + v (65n_ — 102n + — 148) + v 2 (—36n+ — 72)) 
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(A4) 


Appendix B: Explicit results for the source multipole moments 


We list in this appendix explicit results for the newly computed SS contributions to the 
source moments, the computation of which is described in Section II V B1 We recall that the 
brackets indicate the STF projection. 

For the mass quadrupole moment, we obtain 
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(Bla) 


with 

it = - S Ki S j> (5k_ + «+) + 4 S Ki E j> vK- + E <l E 3>u(5k- - k+) , 
ig 0 = 29 S <i S ,J ' V (5k_ - k+) + 58 S Ki E j> v 2 ((«_ - 5^ + ) - 2 vkJ) 

+ E<*E' 7> u 2 (29(5k_ — k + ) + u{— 295k_ + 87 n + + 140)) + 66 (. Sv)S <l v : ’ > (k + — 5n_) 
+ 66 (E v)S <l v^ > ((5k+ — K-) + 2vnJ) + 66 (Sv) E < V 7> ((5k + — «_) + 2uk-) 
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+ 6 (Ev)S<V > (11 (k+ - 5k.) + u(115k- - 33k + - 28)) + 22 S 2 v Ki v j> (5/c_ - k+) 

+ 44 (S'E)u < V ?> ((/c_ — 5 k + ) — 2uk-) 

+ 2 E 2 v <l v‘ > (11(5k_ — ac + ) + z/(— ll5/c_ + 33k + + 14)) , 

^6.i = 6 (raS') 2 n < *n J> ((75 k_ + 18k + + 36) — 40z/(k + + 2)) 

+ 2 S' 2 n <ti n? > ((— 115k_ + 8k + — 48) + 30u(k + + 2)) 

+ 6 (nS)(nYj)n <l n^ > ((115k+ + 365 — 11ac_) + Au(— 105k + — 205 + 3/c_)) 

+ 3 (nE) 2 n < *77,- ?> (ll(/c + — 5kJ) + 2i/(135/c_ — 24k + — 36) + 80u 2 (k + + 2)) 

+ 2 (5'E)n < *rP > ((195 k + — 485 — 19k_) + 2^(15 5k + + 305 + 7k_)) 

+ E 2 n < *n- 7> (19 (k+ — 5k_) — 2i/(45k_ + 15k + + 8) — 60u 2 (k + + 2)) 

+ 12 (nS^S^ ((25k_ - 13k+ - 22) + 5z/(k+ + 2)) 

+ 2 (nE)n < *5' J> (5(—95 k + — 25 + 9k_) + 3^(55k + + 105 — 13k_)) 

+ 2 S <l S f> (17 5 k. + 109k + ) 

+ 2 (nS)n < ’ l TP > ((—455k + — 1225 + 45k_) + 3z/(55 k + + 105 — 13/c_)) 

+ 6 (nE)n < *E- J> (15(5 k_ — k + ) + v(—Q5k- + 39k + + 44) — 10^ 2 (k+ + 2)) 

+ 8 5 ,<i E J> (23(5 k+ - k_) - 17^k_) 

+ 2 E <i E J> (46 (k+ - 5k.) + i/(-175k_ - 75k+ + 56)) . (Bib) 


The current quadrupole moment reads at leading order (see also Ref. 
expressions at any multipolar order) 



for leading order 


(■ J ij ) ss = - 2K-S <i e j>ab S a v b + (-3 - 5k _ + k+^V^E^ 

+ (-5k_ + K + )E < V >a6 S' 0 t; f , + (-«_ + 5 k + + 2fc_z/)E < V >a6 E a u 6 ) . (B2) 
Finally, for the mass octupole moment, we find |42j 

(r ik ) ss = ^n <f (-2 K-S’S k> + 2(k + - ■5k_)S 3 ’E'' > + (&+ - k_ + 2n_(/)E i E‘ > ) . 

(B3) 


Appendix C: Correspondence between the spin vector and spin tensor variables 


This appendix provides the link between the spin tensor and the conserved-norm spin 
vector variables which we use to present our PN results. We recall that the spin tensor 
variable S b is the spatial part, in harmonic coordinates, of the spin tensor introduced in 
Section 111 A1 and that the spin vector variable has been defined in Section 111 Cl as S l = S L , 
with S/j, given by Eq. (12.14[i and i being a spatial index referring to the tetrad e a M constructed 
in the same section. 

We display below the SS contributions to the expression of the spin vector in terms of the 
spin tensor, in the general frame. These contributions complete those computed at the SO 
order in Ref. (Hj, Eqs. (B.l) (notice that the spin tensor components there were denoted as 
S lj instead of S lj ). We have 
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~n a l2 S™ (vfS?£ M ) - Sfe iab (n° V2 v b Sf) + -n a 12 S l 2 a (v“S?£ abi ) 


+S?e iab (n a l2 v b 2 S?) 


+ 0(7 ), 


(Cl) 


where we have authorized the repetition of indices appearing in scalar quantities enclosed 
with parenthesis. At this order, there appear S 1 S 2 terms only and thus no Sf, S 2 terms. 


Appendix D: Equivalence with ADM results for the dynamics 


In this 
in the AI 
of ADM { 

restrict ourselves to the comparison of our findings with the ADM ones, in line with our 
previous works. 

The two results have been obtained in different gauges and the spin variables differ in 
their definition. It is thus important to take properly into account the transformation of the 
particle positions and spins from one formalism to the other. In the following, we will denote 
the ADM variables with an overbar and resort to the convenient notation If a — Pa/ 171 a■ Let 
us now introduce the contact transformation Ya(x,p, S ) and the rotation vector 6 A (x,p, S) 
such that the harmonic variables are related to the ADM ones by 

y A = Y A (x,p,S) + 0( 7), (Dla) 

S A = S A + O a (x,p,S) xS a + 0(6). (Dlb) 


appendix, we compare our results for the dynamics with those previously obtained 

As the equivalence 
we actually 


62|, 170H73] and EFT [21, 


38, 39, 74 


approaches. 

md EFT results has been already demonstrated in Refs. j36|, 


The ADM spin variables and ours have the same Euclidean norm Sa • Sa — Sa • Sa — s\, 
which is precisely the conserved norm introduced in Section III AI Since the first corrections 
enter as 6a = 0(4), we see that the transformation for the spins necessarily takes this form. 

Now, if we denote by Aa(x,p, S) and Qa(x,p, S) the function that converts to ADM 
variables the harmonic-coordinate acceleration and precession vector, and by Qa the pre¬ 
cession vector of the ADM spins, such that dS^/df = ft a x Sa, the two relations to impose 
for the dynamics to be equivalent are 


Aa — {{TA, 77adm} , TLadm} + 0(7 ), 

{6a, Hadm} + 6A X n A = (x, p, S) - n A + 0(6 ), 


(D2a) 

(D2b) 


where T/adm is the ADM Hamiltonian (which can be found for instance in Section 6.2 of 


Ref. 62]) and {,} is the usual Poisson brackets extended to spin variables. Here the term 


6 A x ft a is actually negligible, for 6 4 = 0(2) and ft a = 0(4). 

We find that there are no contributions at leading order in the transformations (IDID . i.e. 
(L^ss — 0(6) and (0a)so = 0(5). Using the method of undetermined coefficients then 
leads to a unique solution for the higher-order terms in the transformations. For the rotation 
vector we obtain 
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+ 0(7), 


(D3) 


with (adopting the same notations as in the rest of the paper for scalar products) 
6/ 1 = ——^-7Ti2 [(T^Si) + (ni27f 2 )(ni 2 S'i)] - ^7f 2 (n r2 5' l ) - ^Tfi(ni 2 Si) 
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^l’ 2 — — 2 ) + 2 ni2 [( 7r 2‘^ 2 ) ~~ 3( n 12 7r 2)(^12'S , 2)] + 7I ’ 2 ( ri 12*S’2) • (D4) 

We recall that SO terms in 9 actually correspond to SS effects in the dynamics. For the 
contact transformation, we arrive at the simple expression 

771 

(U) ss = 2^k [* l( " 12Sl) “ + °< 8 )' < D5 ) 


The relevant NS and SO contributions to these transformations are given for instance in 


Ref. 24] and Refs. 24, [25 


The existence of a solution relating our variables to the ADM ones validates our results, 
the problem of finding such transformations being largely over-constrained. 
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